Cluster Growth in two- and three-dimensional Granular Gases 
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Dissipation in granular media leads to interesting phenomena such as cluster formation and crys- 
tallization in non-equilibrium dynamical states. The freely cooling system is examined concerning 
the energy decay and the cluster evolution in time, both in two and three dimensions. We also 
suggest an interpretation of the 3D cluster growth in terms of percolation theory, but this point 
deserves further study. 
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I. INTRODUCTION 



Granular media are interestingmulti-particle systems 
with a rich phenomenology Q, 0,13, 0, . They can form 
a hybrid state between a fluid and a solid, where the be- 
havior is controlled by the balance between energy input 
and energy dissipation. Energy input leads to a reduc- 
tion of the density due to more collisions and increasing 
pressure, so that the material can flow. In the absence 
of energy input, e.g., in freely cooling systems, granu- 
lar materials become denser, i. e. they locally "solidify" 
due to dissipation. Because of mass conservation, the lo- 
cal densification is accompanied by a density decrease in 
other parts of the system, giving rise to complex patterns 
and structures, with an interesting time evolution. How- 
ever, theoretical approaches are non-classical and appear 
often extremely difficult, so there is still active research 
directed towards the better understanding of granular 
media. 

The subject of this paper is the pattern formation via 
clustering in a dissipative, freely cooling granular gas 
H 0, B E 113 IH 111- The basic idea of clustering is 
that in an initially homogeneous freely cooling granular 
gas, fluctuations in density, velocity, and temperature 
cause a position dependent energy loss. Due to strong 
locally inhomogeneous dissipation, pressure and energy 
drop rapidly and material moves from "hot" to "cold" 
regions, leading to even stronger dissipation and thus 
causing the density instability with ever growing clusters 
until eventually, clusters reach system size (see Fig.^l. 

We investigate this clustering instability with respect 
to different dissipation rates and different system sizes. 
Even though many of our findings are empirical, we at- 
tempt to reduce the complexity of the evolution of the 
system by use of simple scaling laws. 

In section [H] we explain in detail the simulation ap- 
proach. The results of our numerical experiments are 
discussed in section IIIII Finally, in section IIVI a short 
summary and a discussion are given. 




FIG. 1: Density distribution in a snapshot of a 3D system 
with 512000 particles, volume fraction p — 0.25, and restitu- 
tion coefficient r = 0.3. 

Each sphere represents a cell with about 40 particles in av- 
erage. The size of the spheres is proportional to the local 
density; very small spheres (corresponding to low densities) 
are omitted. 



II. SIMULATION DETAILS 

A granular gas can be idealized as an ensemble of hard 
spheres in which the energy loss that accompanies the 
collision of macroscopic particles is modeled with a sin- 
gle coefficient of restitution r. In the simplest case the 
particles are identical in size and mass and there are no 
inter-particle forces between collisions. 

Details about initial and boundary conditions are given 
m section llTAl The microscopic dynamics of the mo- 
tion and the collision of the particles is discussed in sec- 
tion ^TO] and the simulation method is explained in sec- 
tion lll Gl Section Hi Dl deals with the inelastic collapse, a 
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problematic artefact of the hard sphere model with dis- 
sipation. 



A. Initial and Boundary Conditions 

The simulation volume consists of a box with equal 
side length and periodic boundary conditions in two or 
three dimensions. 

An initial state with random particle positions and ve- 
locities is prepared in the following way: The particles 
first sit on a regular lattice and have a Maxwellian ve- 
locity distribution with a total momentum of zero. Then 
the simulation is started without dissipation and runs for 
about 10 2 collisions per particle. This state is now used 
as initial configuration for the dissipative simulations. 



B. Microscopic Dynamics 

Between collisions no forces act upon the particles and 
they move at constant velocity. 

The particles are idealized as hard spheres. This means 
that collisions take infinitesimal time and involve only 
two particles. Conservation of momentum leads to the 
collision rule [EE! 

1 + T / ~ \ ~ 

v 'i/2 = v i/2 T — ^ — [k • (vi — v 2 )J k . (1) 

where a prime indicates the velocities v after the collision, 
k is a unit vector pointing along the line of centers from 
particle 1 to particle 2, and r is the coefficient of resti- 
tution. The relative tangential velocity does not change 
during a collision, the relative normal velocity changes 
its sign and is reduced by a factor 1 — r. Energy dissi- 
pation is proportional to A = 1 — r 2 , so that the elastic 
limit r = 1 implies A = 0, i. e. no dissipation, while r < 1 
implies A > 0. 



The algorithm processes the events one after another. 
After a collision the positions and velocities of the two 
involved particles are updated; the state of all other par- 
ticles remains unchanged. For the two colliding particles, 
new events are calculated and the next future event is 
stored in the event priority queue for both particles. The 
next event is obtained from the priority queue, the new 
positions and velocities after the collision for the collision 
partners are updated, and so on. Neighborhood search is 
enhanced with standard linked cell methods where 
the cell change of a particle is treated as a new event t ype . 
The details of the algorithm can be found in [Til Il5lll7| . 



D. Avoiding the Inelastic Collapse with the TC 
Model 

Our model makes use of hard spheres with an infinitely 
stiff interaction potential. This is an idealization of real 
physical particles and can lead to the dramatic conse- 
quence of inelastic collapse: An infinite number of colli- 
sions occurs in finite time. This singularity is unphysical, 
of course, and a major drawback for numerical simula- 
tions, too. But it has been shown that one can circum- 
vent this artefact of the model in the following way [18j : 
If two consecutive collisions of a particle happen within 
a small time t c , dissipation is switched off for the second 
collision. This timestep t c corresponds to the duration of 
the contact of physical particles. 

There exist other deterministic and random models 
which prevent inelastic collapse. Even though many of 
them lack a solid theoretical background and physical 
motivation, their details should be insignificant for the 
physical evolution of the system anyway, since only a 
small negligible fraction of the particles in the system is 
involved in the inelastic collapse. For an extensive dis- 
cussion see |Igj . 



C. Event-Driven Molecular Dynamics 

The simulation of hard spheres can be handled effi- 
ciently with event-driven molecular dynamics [T3. IT^ |. 
The collisions are the events which have to be treated 
by the algorithm. Between these collisions the particles 
move on trivial trajectories and so the algorithm can eas- 
ily compute the point of time t\i of the next collision of 
two particles 1 and 2 as 

tl2 = t - ri2 • Vi 2 /U 2 2 

+ vWvi 2 ) 2 -(r 2 2 -d 2 ) W 2 2 /„2 2j (2) 

where v i2 = v 2 (*o) - Vi(t ) and r i2 = r 2 (t ) - ri(i ) 
are the relative velocities and positions of the particles 
at time to, and d is the diameter of a particle. 



III. NUMERICAL EXPERIMENTS 

The simulation is started from a homogeneous system, 
prepared as described in section Til Al Depending on the 
dissipation A, the density p, and the number N of par- 
ticles, the system remains in the homogeneous cooling 
state (HCS) for some time, until clustering starts and 
the system becomes inhomogeneous. 

We discuss the evolution of the kinetic energy and the 
collision frequency of the system in section ITlI Al In sec- 
tion ^^B] we investigate the clustering by means of ap- 
propriate measures of the cluster size distribution. In 
section UlI CI and IIII Dl we encounter interesting parallels 
with percolation theory and discuss several critical expo- 
nents. 
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FIG. 2: Decay of the kinetic energy E (left) and the collision frequency / c (right) plotted against scaled time r in a 2D system 
with N = 316 2 = 99856 particles, volume fraction p — 0.25, and different restitution coefficients r (increasing r from top to 
bottom). The thick solid lines correspond to the theoretical predictions as given in the inset. 
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FIG. 3: Decay of the kinetic energy E (left) and the collision frequency f c (right) plotted against scaled time r in a 3D system 
with N = 80 3 = 512000 particles, volume fraction p — 0.25, and different restitution coefficients r (increasing r from top to 
bottom). The thick solid lines correspond to the theoretical predictions as given in the inset. 



A. Kinetic Energy and Collision Frequency 

Dissipativc collisions lead to a decay of the kinetic en- 
ergy and the collision frequency (see Fig. |2 and Fig.|3J). 
From these figures three different regimes can be clearly 
distinguished. They are: 

(1) The homogenous cooling state (HCS) at the be- 
ginning, when no clusters have formed yet, is well un- 
derstood [T^ . |20| . The decay of the kinetic energy E is 
governed by the equation 



E(r) = 



E(0) 



(1 



(3) 



with the scaled time r — ^-r-- D is the dimension 
of the system and tE is the initial Enskog collision time 
t e = {\fTrd) J ' {2\/15vpgd{p))- v is the mean velocity of a 



particle, d its diameter, p is the volume fraction, and g c i 
is the contact probability. In 2D gd{p) = (1 — 7p/16)/(l — 
pf and in 3D > g d (p) - (1 - p/2)/(l - pf . 

The evolution of the collision frequency per particle 
fc(j) with time is given by 



fc(r)=t E \0), 



' E(r) 
25(0) 



(4) 



and it is the natural time scale controlling the evolution 
of the system in the HCS. 

(2) When clusters start to grow, the decays deviate 
from these laws. In the cluster growth regime, the decay 
of the kinetic energy slows down. Furthermore, the col- 
lision frequency fluctuates a lot and can even increase. 
Note that the collision frequency is no longer a natural 
time scale of the system, since it is mainly determined 
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FIG. 4: Growth of the 1st moment Mi (left) and the 2nd moment M2 (right) of the cluster size distribution plotted against 
scaled time r in a 2D system with N — 99856 particles, volume fracion p = 0.25, and different restitution coefficients r 
(increasing r from top to bottom). 




FIG. 5: Growth of the 1st moment Mi (left) and the 2nd moment M2 (right) of the cluster size distribution plotted against 
scaled time r in a 3D system with N = 512000 particles, volume fracion p — 0.25, and different restitution coefficients r 
(increasing r from top to bottom). 



by cluster-cluster collisions, where it increases strongly. 
The deviation from Eq. Q occurs earlier and is more dra- 
matic for larger dissipation A, i. e. smaller r. However, 
the cluster growth regime is characterized by an energy 
decay around E ~ r _1 , independently of r, cf. [ljj. (In 
3D the exponent might be slightly larger than unity; a 
best fit yields E ~ T - 1:L ± 01 .) J n contrast, the evolution 
of the collision frequency with time depends on r dur- 
ing cluster growth. This regime is more distinct for large 
dissipation A and large system sizes N 1 ' and can c. g. 
barely be seen for r = 0.98 in Fig. [3J 

(3) Finally, when the largest cluster in the population 
has reached system size in the /em saturation regime, 
the cooling resembles the homogeneous cooling state in 
so far that E(r) oc r -2 and / c (t) oc t -1 , even if the latter 
still shows large fluctuations. In Fig. [3 this regime is not 



clearly visible, because the simulation times are not long 
enough for this system. 

More details of the 2D situation were discussed in 
ITsj and references therein. Further studies on the three- 
dimensional systems are in progress. 



B. Cluster Growth 

The energy loss of the particles first leads to a reduced 
separation velocity after collision and eventually to the 
formation of clusters. But the definition of a cluster suf- 
fers from the fact that it takes an infinite number of colli- 
sions for the particles to stay in permanent contact with 
each other [|| . So we use the following (geometrical) def- 
inition: Two particles belong to the same cluster, if their 
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distance is smaller than s = 0.1 particle diameters. The 
choice of s is arbitrary and shifts the results; the qualita- 
tive behavior of the quantities discussed below does not 
depend on s, as long as it is neither too small nor too 
large 0. 

The moments Mk of the cluster size distribution are 
defined as 

M k :=— Vi'n,, (5) 

where n c denotes the total number of clusters and rii the 
number of clusters of size i. 

In many cases there are a lot of small clusters and one 
large cluster of size N x , which contains a macroscopic 
fraction m x :— N x /N of the total number N of particles. 
Therefore, we also define reduced moments ML, which do 
not include the largest cluster. 

In Fig. 0] and Fig. [5] the growth of the clusters can be 
seen on the basis of the first and second moments. Af- 
ter several collisions, particles start to cluster and the 
moments of the cluster size distribution grow until they 
reach their "saturation" value • A numerical anal- 
ysis reveals that the increase in M\ and M2 is mainly 
due to one large cluster which grows until it reaches its 
maximum size. In this final state this cluster contains a 
macroscopic fraction m x of the particles (see Fig. 0). 
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FIG. 6: Growth of the large cluster m x in a 3D system with 
N particles, volume fraction p — 0.25, and a restitution coef- 
ficient r = 0.75. 

The onset time of cluster growth and also the final size 
m x of the large cluster depend strongly on the restitu- 
tion coefficient r (see Fig. 0] and Fig. 0). At low dissi- 
pation rates, for a long time nothing interesting happens 
and finally small and strongly fluctuating clusters ap- 
pear. High dissipation leads to almost immediate cluster 
growth and a very large cluster at last. On the other 
hand, the system size N does not seem to affect the be- 
havior of the system provided that N is not too small 
(see Fig. |SJ . 



It has been shown 0, 0, [2l| , that for infinite system 
size the system is always unstable to the formation of 
clusters, whereas for a finite system size the density p and 
the dissipation A = 1 — r 2 must not be too small. With 
our choice N > 10 5 and p = 0.25 we expect no cluster 
formation below the critical dissipation A c = 10 -4 in our 
2D system and A c = 10~ 2 in our 3D system j8j. This is 
in good agreement with our results for the 3D system, 
where cluster formation still can be seen with dissipation 
A = 0.04, but not for A = 0.004 (see Fig.©. 

As can be seen in Fig.0(left) the onset r c of the growth 
of the large cluster in 2D does not depend on the restitu- 
tion coefficient r explicitly. (Implicitly, r c is dependent 
on r via rescaled time, of course.) In contrast t c does 
strongly depend on r in a 3D system (see Fig. (right)). 

An empirical formula, which gives the dependency 
rather accurately over orders of magnitude in A, is 

T C (A) = T 1 + C (i-^, (6) 

where T\ = 0.24 ± 0.02 and c = 36 ± 5 are fit parameters 
(for s = 0.1). 

In the limit A — > 1 , one has r c — > t± , where t\ is propor- 
tional to the rescaled collision frequency and corresponds 
to t = 2Dt\IeI\ ~ 1.5 tE- This constant stems from the 
fact that the particles need at least one collision to clus- 
ter, however, it also depends on s. In the limit A —> A c 
the time t c diverges: the particles never start to cluster 

M 

of magnitude in Rescaling time in Fig. [7] (right) ac- 
cording to Eq. © results in Fig. [H] (right). There, the on- 
set of the clustering happens approximately at the same 
rescaled time, whereas these times differ by 4 orders of 
magnitude in Fig. (right). The fluctuations in r c are 
not systematically dependent on r. 

Another difference between the 2D and 3D system in 
Fig.0is the fact that the cluster growth in 2D starts very 
smoothly, whereas the beginning of the clustering is very 
sharp in 3D. 

C. Critical exponents in 3D 

As the moments of the cluster size distribution are 
clearly dominated by the large cluster, we will now study 
the reduced second moment M%, which does not include 
the large cluster. Fig. [S] shows that M' 2 is small most 
of the time, except for a peak at r c . This gives us a 
clean definition of r c . But what is more interesting, par- 
allels to percolation theory arise. In percolation theory 
|22| one studies the scaling behavior of certain quanti- 
ties depending on the occupation probability p around 
the percolation threshold p c . One result for the reduced 
second moment M 2 (p) is, e.g., 

Mfo)~\ p - Pc \-r , (7) 

where 7 is a universal critical exponent. In order to trans- 
fer this result to cluster growth, we replace p with r and 
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FIG. 7: Growth of the large cluster m x in 2D (left) and 3D (right) for different restitution coefficients 





FIG. 8: (left) Scaling behavior of the onset r c of the cluster growth in 3D depending on the dissipation rate A. The curve is 
Eq. @. 

(right) With rescaled time r/r c according to Eq. © the curves of the growth of the large cluster in 3D almost collapse. 



p c with r c : 



M' 2 {t) 



(8) 



Indeed we find in Fig. the same power law behaviour. 
All data for different system sizes collapse on the same 
master-curve and even the exponent 7 = 1.8 ± 0.1 is in 
agreement with the 3D percolation problem result 7 = 
1.80 after the beginning of clustering. 

Now, if we study the amplitude of this peak, percola- 
tion theory tells us that the maximum sits at 



Pn 



P. 



cil-aL- 1 /") 



(9) 



where L is the system size, a is a constant, and v is 
another critical exponent. Thus we expect: 



M' 2 { Pmax ) ~ 



(10) 



As the system size L ~ TV 1 / 3 , 
peak in Fig. scales as 



the maximum value of the 



M'{N)~W'* V 



(11) 



Here we can also verify the power law behavior (see 
Fig. 1101) . Our simulations yield the result 7/3^ = 0.77 ± 
0.09, which leads to v = 0.78 ± 0.14. Within the rather 
large margin of errors this result is close to the 3D per- 
colation problem v = 0.88, too. 

A third exponent q [23 is given by the cluster size 
population at time r c : 



n 4 (T c ) 



(12) 



where rij is the number of clusters of size i. Fig. 1111 
yields q = 2.2 ± 0.2. This result is in agreement with 
the prediction of percolation theory q — 2.2, too. 

Is this similarity of clustering in 3D and percolation 
coincidence or is there a deeper reason? 
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FIG. 9: Reduced 2nd moment M'i of the cluster size distribution in 3D systems with different N, volume fraction p = 0.25, 
and a restitution coefficient r — 0.6. The data are averaged over 10-20 different simulation runs, 
(left) The curves show a peak at the onset r c of the cluster growth. 

(right) Scaling behavior of M^(r) against \t — T c \; each data curve has two branches, solid and dotted lines connect data points 
for r > t c and r < r c , respectively. The additional line has a slope of -1.8. 




FIG. 10: Maximal amplitude of the reduced 2nd moment M2 
of the cluster size distribution (see Fig. |§J. The line has a 
slope of 0.77. 



FIG. 11: The points give the numbers n; of clusters of size i 
at time r c in a 3D simulation with N = 512000 particles and 
a restitution coefficient r = 0.75. The line has a slope of -2.2. 



D. Consequences 

Percolation theory makes many universal predictions 
for randomly disordered systems. The particles in a 
granular gas form such a disordered system. There have 
been other attempts to apply percolation theory to gran- 
ular systems j2^| , and also more advanced phase ordering 
models have been used to parallel the clustering dynam- 
ics ^1^2- Our results indicate that such attemps can be 
justified, even though the percolation problem is purely 
static, whereas the clustering of granular gases also in- 
volves the dynamics, i. e. momentum and energy. 

The only thing that might seem strange at first sight is 
the derivation of Eq. JSJ. Why can we replace the purely 



static quantity p with time r? 

In order to introduce the occupation probability p in 
the clustering problem, we define it as 



p(r) 



M 1 {t)-M 1 {t*) 
Mi(t) 



(13) 



where r* < t c is a time shortly after the first few colli- 
sions. Now, we examine p(r) around r c . Fig. 1121 (right) 
shows that in 3D, p is almost proportional to r. If we 
make use of this linear relation and insert it in Eq. J7|), 
we arrive at the postulated formula Eq. (JSJ. 

In contrast, in our 2D system the growth of the large 
cluster happens at a rather large time r c . (Note, that r c 
is not very clearly defined in the 2D case, because the 
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FIG. 12: "Occupation probability" p against scaled time r/r e . 
(left) In 2D the data curves are highly non-linear, 
(right) In 3D the data curves are almost linear. 



onset of the cluster growth is very smooth.) In this late 
stage, the relation between p and r is highly non-linear, 
see Fig.ll2l(lcft). Thus we could not find universal critical 
exponents for r there. However, it might be possible to 
find a linear relation for a different set of parameters, 
e. g. for a different volume fraction, too. Whether a 2D 
system and a 3D system with the same parameters (like, 
e.g., volume fraction) are in fact comparable is another 
open issue. 



IV. SUMMARY AND DISCUSSION 

The evolution of freely cooling granular systems can 
be divided into three regimes. First, the system is in the 
homogeneous cooling state (HCS). Then, in the cluster- 
growth regime, clusters begin to develop and grow. Fi- 
nally, in the saturation regime, the clusters merge to prac- 
tically one large cluster, which grows until it reaches sys- 
tem size. Besides the macroscopic fraction of particles in 
the large cluster, there are still many small clusters with 
interesting statistics. 

In the HCS, the decay of the kinetic energy E and 
the collision frequency f c can be described by the sim- 
ple analytical expressions E(t) ~ (1 + t)~ 2 and / c (t) ~ 
(1 + t) . The collision frequency is the natural time 
scale here, mainly determined by the density and the 
dissipation rate of the system. For strong dissipation, at 
short range, particles already stay closer together after 
only one collision, so that the moments of the size distri- 
bution change rather early due to a (short-range) change 
in the radial pair-distribution (data not shown). 

In the cluster- growth regime, the collision frequency 
shows large fluctuations because of cluster-cluster col- 
lisions and cannot be predicted during cluster growth, 
because it changes erratically and discontinuously. The 
energy decay is characterized by E ~ r , where the ac- 



curacy of the exponent is limited in 3D due to the com- 
paratively short duration of the cluster growth regime. 
However, this regime shows interesting differences be- 
tween two and three dimensions. In 3D, cluster growth 
can be described by a power law behavior with the same 
critical exponents as in percolation theory. The onset 
of this cluster growth t c is very sharp and does strongly 
depend on dissipation A. In contrast, we could not find 
a similar behavior in 2D, where the onset of the cluster 
growth is very smooth and depends on dissipation A only 
implicitly. 

When cluster growth has reached a dynamic equilib- 
rium in the saturation regime, the system is dominated by 
one large cluster which contains a macroscopic fraction 
of the system. Note that this regime still has an interest- 
ing dynamics and smaller clusters interact with the large 
cluster. Neither the small clusters nor the large one are 
static and the latter loses or eats up particles or smaller 
clusters. Kinetic energy and collision frequency still fluc- 
tuate, but are governed by the equations E(t) ~ t~ 2 and 
Ic{t) ~ t -1 . This means the evolution in time is similar 
to the homogeneous cooling state. 

A very interesting observation is the similarity of clus- 
tering in 3D and percolation with respect to the critical 
exponents. Even for the dynamics of cluster growth we 
have found the same exponents as for the occupation 
probability in percolation theory. We have provided a 
tentative explanation which is based on the linear rela- 
tion of the occupation probability p and time r around 
the onset r c of cluster growth in 3D. In 2D we did not 
find a power law behavior because the growth of the large 
cluster happens when the linear relation between p and 
t has disappeared — at least for the given set of param- 
eters. 

However, there still remain some open questions: Why 
does t c depend on r in 3D, but not in 2D? Why does the 
growth of the large cluster proceed smoothly in 2D and 
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rather sharply in 3D? And why are all these differences 
present while the energy decay is proportional to t _1 in 
both cases? Detailed studies of these questions and also 
of the cluster size probability distribution are in progress 
as well as a more systematic study of the cluster definition 
in 3D. 



Acknowledgments 

This research was supported by the DFG, the SFB 382, 
the SFB 404, and LU 450/9-1. We thank Hans Herrmann 
and Sean McNamara for helpful discussions. 



[1] H. J. Herrmann, J.-P. Hovi, and S. Luding, eds., Physics 
of dry granular media - NATO ASI Series E 350 (Kluwer 
Academic Publishers, Dordrecht, 1998). 

[2] T. Poschel and S. Luding, eds., Granular Gases 
(Springer, Berlin, 2001), lecture Notes in Physics 564. 

[3] P. A. Vermeer, S. Diebels, W. Ehlers, H. J. Herrmann, 
S. Luding, and E. Ramm, eds., Continuous and Dis- 
continuous Modelling of Cohesive Frictional Materials 
(Springer, Berlin, 2001), lecture Notes in Physics 568. 

[4] Y. Kishino, ed., Powders & Grains 2001 (Balkema, Rot- 
terdam, 2001). 

[5] R. Cafiero, S. Luding, and H. J. Herrmann, Phys. Rev. 

Lett. 84, 6014 (2000). 
[6] I. Goldhirsch and G. Zanetti, Phys. Rev. Lett. 70, 1619 

(1993). 

[7] I. Goldhirsch, M.-L. Tan, and G. Zanetti, Journal of Sci- 
entific Computing 8, 1 (1993). 

[8] S. McNamara and W. R. Young, Phys. Rev. E 53, 5089 
(1996). 

[9] S. Luding and H. J. Herrmann, Chaos 9, 673 (1999). 
[10] X. B. Nie, E. Ben-Nairn, and S. Y. Chen, Phys. Rev. 

Lett. 89 (2002). 
[11] S. K. Das and S. Puri, Physica A 318, 55 (2003). 
[12] S. K. Das and S. Puri, Europhys. Letters 61, 749 (2003). 
[13] S. Luding, in Physics of dry granular media - NATO 

ASI Series E350, edited by H. J. Herrmann, J.-P. Hovi, 

and S. Luding (Kluwer Academic Publishers, Dordrecht, 

1998), p. 285, also see references therein. 



[14] B. D. Lubachevsky, J. Comput. Phys. 94, 255 (1991). 

[15] S. Miller and S. Luding, J. Comp. Phys. 193, 306 (2004), 
physics/0302002. 

[16] M. P. Allen and D. J. Tildesley, Computer Simulation of 
Liquids (Oxford University Press, Oxford, 1987). 

[17] B. D. Lubachevsky, Int. J. Comput. Simul. 2, 372 (1992). 

[18] S. Luding and S. McNamara, Granular Matter 1, 113 
(1998), cond-mat/9810009. 

[19] S. Luding, M. Huthmann, S. McNamara, and A. Zip- 
pelius, Phys. Rev. E 58, 3416 (1998). 

[20] P. K. Haff, J. Fluid Mech. 134, 401 (1983). 

[21] N. Sela and I. Goldhirsch, Phys. Fluids 7, 507 (1995). 

[22] D. Stauffer and A. Aharony, Introduction to Percolation 
Theory (Taylor & Francis, London, 1994), 2nd ed. 

[23] J. Tobochnik, Phys. Rev. E 60, 7137 (1999). 

[24] For a detailed study of different s in 2D see [IJ; in 3D, 
the s dependence of some quantities seems to be stronger 
than for others. 

[25] Saturation means that there is a transition to a regime 
with a much slower change of the moments. The dynam- 
ics is not finished, however, since the large cluster still 
can collect particles or break into pieces. 

[26] Since we focused on strong dissipation in this study, 
rather than the limit A — > A c , we can not draw a con- 
clusion on the functional behavior from our limited data. 

[27] The common notation for this exponent is r. In order to 
avoid name conflicts we call it q. 



